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ABSTRACT 


Two methods for calculating linear frequency domain unsteady aerodynamic 
coefficients from a time-marching full -potential cascade solver are developed and 
verified. In the first method, the Influence Coefficient method, solutions to 
elemental problems are superposed to obtain the solutions for a cascade in which 
all blades are vibrating with a constant interblade phase angle. The elemental 
problem consists of a single blade in the cascade oscillating while the other blades 
remain stationary. In the second method, the Pulse Response method, the 
response to the transient motion of a blade is used to calculate influence 
coefficients. This is done by calculating the Fourier transforms of the blade 
motion and the response. Both methods are validated by comparison with the 
Harmonic Oscillation method, in which all the airfoils are oscillated, and are 
found to give accurate results. The aerodynamic coefficients obtained from these 
methods are used for frequency domain flutter calculations involving a typical 
section blade structural model. Flutter calculations are performed for two 
examples over a range of subsonic Mach numbers using both flat plates and 
actual airfoils. 
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sonic velocity 

elastic axis position from midchord, see Figure la. 
airfoil semi-chord 
airfoil chord 
lift coefficient 

moment coefficient about elastic axis 
blade motion 
response to fit) 

Fourier transform of fit) 

Fourier transform of Fit) 

response to unit step function in blade displacement 

Fourier transform of F&it) 

cascade gap or spacing, see Figure lb. 

gap-to-chord ratio 

plunging displacement, normal to airfoil chord 
imaginary unit, V-T 
unit matrix 

moment of inertia about elastic axis 
imaginary part of { } 
reduced frequency, k c =icc/Mooaoo 
spring constant for plunging 
spring constant for pitching 
tali i la a 

frequency domain unsteady aerodynamic coefficients 
mass of typical section 
Mach number 

number of blades in the cascade 

matrix defined in Equation (6) 

complex representation of force on 0 th blade 

influence coefficient, force on blade due to oscillation of j th blade 

lift and moment about elastic axis 

interblade phase angle index 

radius of gyration about elastic axis, (7 a /m6 2 ) 1/2 

real part of { } 

blade index 

static unbalance 

time 

fluid velocity in x, y directions 
Cartesian coordinates 

distance between elastic axis and center of mass, x a - S a hiib 
eigenvector in Equation (6) 
reduced velocity, V*= Moo Oo°/bco a 
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pitching displacement about elastic axis 

ratio of specific heats 

unit step function 

stagger angle 

eigenvalue in Equation (6) 

mass ratio, p ■ m/npoob 2 

real part of the eigenvalue (damping ratio) in Equation (7) 

imaginary part of the eigenvalue (damped frequency) in Equation (7) 

transformed spatial coordinates 

fluid (air) density 

interblade phase angle 

nondimensional time, r= a oo tic 

velocity potential 

oscillation frequency 

uncoupled natural frequency for bending (plunging), (Oh “ (Kh,lm ) l/2 
uncoupled natural frequency for torsion (pitching), u) a “ ( K a ll a ) l/2 

d( )ldt 

( ) at far upstream conditions 
( ) at flutter 

amplitude of harmonic variable ( ) 

( ) in the i 4h interblade phase angle mode 
( ) for the s th blade 



INTRODUCTION 


Accurate prediction of flutter in turbomachinery components has been a 
primary research goal of aeroelasticians. The development of advanced turboprop 
(propfan) engines has led to a renewed interest in the study of flutter in bladed 
disks. Propfans, which have been designed to operate at high subsonic cruise 
speeds, have encountered flutter in the early part of their development . 1 

Two fundamental approaches have been used in flutter calculations — 
frequency domain analysis and time domain analysis. The frequency domain 
method is the conventional approach in turbomachinery aeroelasticity. Recently, 
the time domain method, which has found wide application in aircraft or fixed- 
wing aeroelasticity, has also been applied to bladed disk or cascade flutter 
calculations . 2 - 8 The frequency domain approach has been generally preferred 
over the time domain approach because the frequency domain analysis requires 
only unsteady aerodynamic forces corresponding to harmonic blade motion. The 
time domain method, on the other hand, requires aerodynamic forces 
corresponding to arbitrary blade motions. The aerodynamic forces corresponding 
to harmonic blade motions are considerably easier to determine than those 
corresponding to arbitrary blade motions because the explicit time -dependence 
can be eliminated from the governing equations, thus reducing the number of 
independent variables. Further, in a linear flutter analysis, in which only the 
onset of flutter is to be predicted, blade motions can be assumed to be 
infinitesimally small. This allows considerable simplification in the aerodynamic 
analysis since it allows perturbation techniques to be used. For linear flutter 
analysis, the frequency domain and time domain methods yield the same 
results . 2 

The unsteady aerodynamic analyses that have been used in flutter 
calculations have, for the most part, been restricted to inviscid flow theories. 
Almost all the aerodynamic analyses used in flutter prediction have been based on 
linearized potential theories. The earliest analyses 4 - 6 assumed that the unsteady 
component of the flow was a small harmonic perturbation about a uniform, 
steady mean flow; this is equivalent to neglecting the blade thickness and camber. 
Over the past decade, unsteady linearized analyses have been developed 6 ' 8 that 
include the effects of non-uniform mean flow due to blade shape. These analyses 
split the flowfield into steady and unsteady components. The assumption that the 
unsteady flow component is a small harmonic perturbation about the steady flow 
allows the elimination of the time variable from the governing equations. The 
unsteady component of the flow is described by a linear equation with coefficients 
that depend on the mean steady flowfield. 

Recently, unsteady aerodynamic analyses based on the time domain solution 
of the full-potential and Euler equations have been developed . 9 - 10 These methods 
can calculate the unsteady flowfield in a cascade for blade motions that are 
arbitrary functions of time. The blade motions can be prescribed or determined 
from a structural model. Accordingly, these flow solvers have been coupled with 
structural models for time domain flutter analysis . 2 - 11 Also, by specifying the 
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blade motions to be simple harmonic (in time), it is possible 12 to calculate the 
aerodynamic forces required in a frequency domain flutter analysis. In this case, 
although the blade motion is prescribed to be harmonic, no assumption is made 
regarding the linearity of the aerodynamic response (lift and moment). By 
selecting an amplitude that is sufficiently small to ensure a linear response, it is 
possible to obtain the same information that was obtained by linearizing the 
original unsteady equations (by assuming infinitesimal harmonic motions). But, 
this is not a computationally efficient method for obtaining this type of 
aerodynamic data for the following reasons. Firstly, the additional independent 
variable, namely time, adds to the computational cost. Secondly, multiple 
interblade passages must be included in the calculation to account for the phase 
lag between the motion of adjacent blades and this increases computational effort 
considerably. In the present work, alternate methods are developed to calculate 
the frequency domain unsteady aerodynamic coefficients in a more efficient 
manner. 

An approach has been developed in the present study that allows the time 
domain full-potential flow solver 9 to be used efficiently to determine aerodynamic 
data corresponding to small amplitude harmonic blade motions for use in 
frequency domain flutter calculations. The conventional approach to calculating 
harmonic aerodynamic data from a tune domain analysis is to specify the blade 
motions to be simple harmonic and to decompose the resulting time-dependent 
response into Fourier components. This method, in which all the blades are 
oscillated with a specified frequency and and interblade phase angle, is referred to 
as the Harmonic Oscillation method. The Harmonic Oscillation method does not 
exploit the linearity of the unsteady flow problem for small amplitudes of motion. 
In the present work, the principle of superposition is used in the Influence 
Coefficient method to determine the aerodynamic forces for different values of 
interblade phase angle by summing the solutions to elemental problems. The 
elemental problem consists of a cascade with one blade oscillating in harmonic 
motion while the other blades remain stationary. Next, the Pulse Response 
method is combined with the Influence Coefficient method to determine the 
aerodynamic coefficients for harmonic blade motions at different values of 
frequency and different values of interblade phase angle. This is done by giving 
one of the blades in the cascade a transient motion. The resulting transient forces 
on all the blades are Fourier transformed and combined to determine the 
aerodynamic forces at a given frequency and interblade phase angle. These 
methods are validated by comparison with the Harmonic Oscillation method. 
Finally, the aerodynamic data thus obtained is used for flutter calculations. 
Flutter calculations are done with a typical section structural model (lumped 
parameter representation) for each blade. Two examples are considered, one with 
a single-degree-of-freedom typical section and the other with a two-degrees- of- 
freedom typical section. The structural parameters and airfoil shape used in the 
second example are representative of the SR5 propfan. 13 
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ANALYSIS 


Flutter analysis of a fluid -structure system requires an aeroelastic 
representation of the system. Accordingly, the structural model, the aerodynamic 
model and the aeroelastic stability analysis are described in this section. 


Structural model 

Each blade of the rotor (bladed disk) is represented by a typical section 
model 14 . The dynamic behavior of the blade is represented by a rigid two- 
dimensional airfoil section taken from the three-quarter span location of the 
blade. The typical section is allowed two degrees of freedom — plunging and 
pitching as shown in Figure la. The plunging motion is in a direction 
perpendicular to the airfoil chord and models the bending of the blade; the 
pitching motion about the elastic axis simulates the torsion of the blade. The 
equations of motion for a single blade are: 


m h 

+ S a a 

+ K h h 

= Q h 

(1) 

S a h 

+ I a a 

+ K a a 

- Qa 

(2) 


where h is the plunging (bending) displacement, a is the pitching (torsion) 
displacement, m is the airfoil mass, la is the moment of inertia, S a is the static 
unbalance, K/ t and K a are the spring constants for plunging and pitching, 
respectively; Qh and Q a are the aerodynamic loads (lift and moment, respectively) 
and the dots over the displacements indicate differentiation with respect to time. 


Aerodynamic model 

The aerodynamic model is based on the unsteady, two-dimensional, full- 
potential equation. The governing equation for irrotational, isentropic flow written 
in conservative form is: 
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In these equations, <p is the velocity potential; p°° , a » and are the density, sonic 

velocity and Mach number, respectively, all evaluated at far upstream conditions. 

The governing equation is transformed to the computational plane where it is 
discretized and solved using a finite-volume scheme. The method of solution 9 * 15 is 
a time-marching scheme that uses approximate factorization at each time level 
with quasi-Newton iterations to maintain time accuracy. The transformation 
from physical coordinates ( x , y, t) to computational coordinates (cf, rj, t) is defined 
at discrete grid points. An H-type grid, shown schematically in Figure lb, is used. 
It consists of £=constant lines that are parallel to the y-axis and constant lines 
that are obtained by interpolating between the upper and lower boundaries of each 
interblade passage. For unsteady calculations, the airfoil positions change with 
time and a new grid is generated at each time level by interpolating between the 
passage boundaries determined from the instantaneous positions of the airfoils. 

The governing equation iB solved as an initial -boundary value problem. For 
steady calculations, a uniform flowfield is used as the initial condition; for 
unsteady calculations, the steady flowfield is used bb the initial condition. The 
boundary conditions are described below. 

Airfoil surface: 

The airfoil surfaces sure treated as impermeable, i.e., the normal velocity of 
the fluid relative to the airfoil surface is zero. Flow tangency is imposed at the 
instantaneous position of the airfoil at each tune step. In order to affect this, a 
new grid which conforms to the airfoil surfaces is generated at each time step. 
Blade motions are thus fully accounted for using moving grids. 

Wake: 

The vorticity that is shed from the trailing edge of the airfoil is assumed to 
remain confined within an infinitely thin region of the flowfield denoted as a wake 
sheet. Wakes are free material surfaces with no pressure difference across them. 
In the present analysis, wake positions are prescribed in advance; continuity of 
pressure and normal velocity is then enforced across the wake. This procedure is 
followed to avoid difficulties related to the tracking of the wake location at each 
time step. The wakes are prescribed to be straight lines extending from the 
trailing edge of each airfoil to the exit boundary. The slope of these lines is selected 
to be the mean of the slopes of the upper and lower surfaces of the airfoil at the 
trailing edge. 

Inlet/ Exit boundaries: 

The inlet and exit boundaries are artificial boundaries that have been 
introduced because the numerical computations cannot be extended to infinite 
distances in the upstream and downstream directions and must be restricted to a 
small, finite region in space. Since these are non-physical boundaries, the 
boundary conditions that are used must allow acoustic waves to pass through 
without reflecting them back into the computational region. Even though it would 
be convenient to prescribe the far upstream/downstream flow velocity or Mach 
number at these boundaries, this is not done because it leads to non-physical wave 
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reflections. In the present analysis, characteristic boundary conditions 16 are used 
at the inlet and exit boundaries. At the inlet/exit boundary, the Riemann 
invariant corresponding to the incoming characteristic (positive characteristic 
with respect to the inward normal) is prescribed based on the far 
upstream/downstream conditions. This allows planar acoustic waves, which are 
at normal incidence to the boundary, to pass through without reflection. However, 
the acoustic waves are generally not planar; also, they are generally not at 
normal incidence to the boundary. Therefore, the use of one-dimensional 
Riemann - invariant boundary conditions results in some wave reflections back 
into the calculation domain. 

Periodic boundaries: 

Periodic conditions are imposed on grid lines extending from the leading edge 
of the airfoil to the inlet boundary; see Figure lb. The periodic conditions are used 
to simulate the fundamental periodicity present on the bladed disk in the 
circumferential direction. 

Additional details concerning the aerodynamic model and the boundary 
conditions can be found in References 9 and 15. 


Aeroelastic Stability 

For a tuned cascade, in which all the blades are identical, the aeroelastic 
modes consist of each of the N individual blades vibrating with equal amplitudes 
with a fixed interblade phase angle between adjacent blades 17 . For the s th blade 
vibrating in the r th interblade phase angle mode, this can be written as: 


,V*l x \^o^\ e i{wt*a r 8) 

\ «s [ 1 «o j 


s - 0, 1, 2,... , AM (4) 


The phase angle between adjacent blades is given as 

o r = 2m' IN r m 0,1,2 AM (5) 

The corresponding aerodynamic forces can be written as linear functions of the 
displacements using the complex-valued frequency domain unsteady 
aerodynamic coefficients Ihh , lha . lah and /«« : 

I Q/isl m b \ - to 2 j ijilir ^0 lb + ih a r a o \ e i (tvt+ a r s) 

\Q as l m ^ 1 \ ^ \ + laar a o j 

where y = mlnp^b 2 is the mass ratio of the blade. The eigenvalue problem for the 
tuned system can be written 18 as: 
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( 6 ) 


where 


[P] = 


M + IjUir 
V (UlJu a ) 

f**a + hhr 


£ *g + jftgr 
£ *j + *ggr 


In the above, x a = S a lmb is the distance between the elastic axis and center of 
mass and r a - da /mb 2 ) 112 is the radius of gyration about the elastic axis, both in 

semi-chord units; a% - C Kh/m) m and u) a m (K a d a ) m are the uncoupled natural 
frequencies for bending and torsion, respectively; note, the subscript ’s' which 
identifies the blade has been dropped. 

For each interblade phase angle mode, the solution of the above eigenvalue 
problem results in two complex eigenvalues of the form 

i JiL = _L_ = /J + iv 

(I (7) 


The real part of the eigenvalue ( /u ) represents the damping ratio and the 

imaginary part ( v ) represents the damped frequency; flutter occurs if /u 1 0 for 
any of the eigenvalues. The stability of each phase angle mode is examined 
separately. The interblade phase angle is fixed at one of the values given by 
Equation (5) and the 2*2 eigenvalue problem is solved. 

However, before the eigenvalue problem can be solved, the aerodynamic 
coefficients must be calculated. For fixed cascade geometric parameters, the 
aerodynamic coefficients are functions of inlet Mach number (Moo), reduced 

frequency of blade vibration (k c ) and interblade phase angle (a r ). The following 
procedure has been adopted for flutter calculations. A value of inlet Mach number 
is selected. A value of reduced frequency is assumed and the aerodynamic 
coefficients Ihh . ha , l ah and l a a are calculated for all values of phase angle in 
Equation (5). The eigenvalue problem is solved for each value of interblade phase 
angle. The reduced frequency is varied until the real part of one of the eigenvalues 

(/ii) becomes zero while the real part of the other eigenvalue is negative. The 
assumed flutter reduced frequency ( k c f ) and the calculated flutter frequency ( v f ) 
are both based on cof. Thus, these two can be combined to eliminate cof and the 

flutter reduced velocity ( V * f) is obtained, viz., V*f = 2 v f /k c f • The critical phase 
angle is identified as the one which results in the lowest flutter velocity. The 
stability of a single-degree-of-freedom pitching system can be inferred entirely 
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from the imaginary part of the moment coefficient, Im{/« a }; flutter occurs if 
Im{/«a}^ 0. 


Aerodynamic Coefficients 

The calculation of the frequency domain unsteady aerodynamic coefficients 
Uhh , ha , lah and / aa ) is described next’ For a given cascade geometry and inlet 
Mach number, these coefficients are required for plunging and pitching motions 
of specified frequency and specified interblade phase angle (restricted by Equation 
(5)). The full-potential solver is based on a time-marching algorithm and three 
different methods of calculating the aerodynamic coefficients are described. 

1. Harmonic Oscillation Method 

In the Harmonic Oscillation (HO) method^, all the N airfoils are oscillated 
with the specified frequency and specified interblade phase angle. The motion of 
the s th airfoil in the cascade is prescribed to be of the form: 

h = h Q sin(aj t+s<J r ) = h 0 Bin(k c M 00 T +s<r r ) for plunging 

or 

a “ a Q sin(t<; t+sa r ) = a 0 sin(fc e r + s <r r ) for pitching 

Calculations are started from a steady flowfield and continued until starting, non- 
periodic, transients have decayed and the flowfield has become periodic in time. 
The lift and moment acting on the reference airfoil (s=0) are calculated at every 
time step and later decomposed into Fourier components. If the amplitude of 
oscillations is sufficiently small, it can be shown by a perturbation analysis that 
the flowfield will have the same harmonic time dependence as the motion 7 . 
Therefore, the lift and moment coefficients (c/ and c m ) on the zeroth blade can be 

represented in terms of a complex quantity Qoas [ Im{Qo} cos(a? t) + 
Re{Q 0 } sin(a) t) ]. The frequency domain unsteady aerodynamic coefficients dhh , 
ha , lah and l aa ) are then obtained by dividing Qo by k c 2 , the amplitude of motion 
and other constants. The harmonic oscillation method requires that calculations 
for each of the AT values of interblade phase angle be done separately. It is possible 
to reduce this computational effort by using the Influence Coefficient method 
described next. 

2. Influence Coefficient Method 

The Influence Coefficient (IC) method is based on the principle of linear 
superposition. Briefly, the solution to a problem is obtained by superposing the 
solutions to the individual elemental problems that comprise the original 
problem. Since the method is based on the principle of linear superposition, it is 
valid only for linear problems. It can be shown that the unsteady part of the 
present problem is linear (governed by a linear differential equation) for 


10 



sufficiently small amplitude of oscillation. 7 It should be emphasized that only the 
unsteady part of the problem is linear; the steady flowfield is described by a 
nonlinear equation. 

Since the quantities of significance are the lift and moment coefficients, the 
following discussion will deal only with a general integrated force quantity. 
However, the results obtained can be extended by analogy to pressure distributions 
or the distribution of other flow variables. Also, complex notation is used for 
convenience; in this notation it is implied that only the imaginary part of the 
complex quantity is to be considered. For example, when the motion is specified to 

be of the form h 0 ^ w * , it is understood that the motion is actually /i 0 sin(u; t). 

The problem to be solved consists of a cascade of N blades in which each blade 
oscillates with a motion of the form sin(o? t+s<r r ), where s is the blade index that 
varies from 0 to AM , and a r is the interblade phase angle given by Equation (5). 
This problem is divided into N elemental problems. The k th elemental problem 
consists of the same cascade of N blades in which the blade oscillates with a 
motion of the form sin(w t) while all other blades remain stationary. The original 
problem and all the elemental problems have solutions that are harmonic 
functions of time. 

For the problem in which all blades oscillate with a motion of the form 
e i(u) t+sa r ) f the forces (c; and c m ) on the 0 th blade can be represented as Qo e iw * \ Qo 
is complex valued to allow the force to lead or lag the motion. The forces on the 0^ 
blade in the k th elemental problem can be represented as Qo,fc e lw Qo,k is 
sometimes referred to as an influence coefficient. Thus, using superposition, the 
following relation can be obtained. 

Qo*'"'- I «o,» (8) 

k-0 

Now, due to the periodicity of the cascade, only the relative positions of the 
oscillating blade and the reference (zeroth) blade are important. That is, the forces 
generated on the 0 th blade due to the oscillation of the k th blade are equal to the 
forces on the blade due to the oscillation of the k+\ th blade, and so on. Thus, 

Qo,k = Q-k , o ■ QN-k , o (9) 

where the periodicity of the cascade of N blades has been invoked again in the last 
step. Thus, the solution to the problem in terms of the influence coefficients can be 
written as: 


N 

Qoe ,w( = £ QN-k, oe ,wt e ,ka r 

k- l 


(10) 


li 


Replacing the influence coefficients Qo,k by the coefficients QN-k , 0 means that 
all the required influence coefficients can be determined simultaneously rather 
than separately. Thus, instead of oscillating the k th blade, calculating the 
pressure history on the zeroth blade and then repeating for all values of k between 
0 and N-l, it is possible to oscillate the zeroth blade and calculate the pressures on 
all the blades simultaneously. This means that the computational effort required 
for the calculation of all the influence coefficients can be reduced by a factor of N 
over the Harmonic Oscillation method. 

It should be noted that the elemental problem used in the present work is 
different from the one used more frequently, in which a single blade in an infinite 
cascade is oscillated 19 * 20 . In the present work, periodic boundary conditions are 
used to simulate an infinite cascade using N blades in the calculations. Thus, 
when a single blade in the iV-blade cascade is oscillated, it corresponds to every 
N tfl blade in the infinite cascade being oscillated. The expressions presented in 
this section are exact for the discrete values of interblade phase angle given by 
Equation (5) and approximate for all other values. Although the summation in 
Equation (8) only extends over a finite number of terms, it does not represent a 
truncation of an infinite sum as long as the value of interblade phase angle is 
restricted by Equation (5). Calculations performed without using periodic 
boundary conditions simulate a finite section of an infinite cascade surrounding 
the blade that is being oscillated. This introduces an approximation into such 
calculations for all values of interblade phase angle; however, this error 
decreases rapidly as the number of blades used in the calculations is increased. 

3. Pulse Response Method 

For a given motion, plunging or pitching, the Harmonic Oscillation method 
and the Influence Coefficient method require separate calculations for each 
oscillation frequency of interest. In order to reduce the computational effort, the 
Pulse Response (PR) method described in this section can be used; this method 
has evolved from the indicial approach that is widely used in many different 
fields. 

Researchers have investigated the indicial approach for aerodynamic 
calculations with isolated airfoils 21 and cascades 22 . Ihe indicial response is the 
response, lift or moment, to a step change in the given mode of motion. From the 
indicial response, the response for any arbitrary motion, specifically harmonic 
motion, can be calculated using Duhamel's superposition integral. Let the time- 
dependence of the blade motion (plunging or pitching) be denoted as fit) and let 
the corresponding response (lift or moment) be denoted as Fit). Let F$ it) denote 

the response corresponding to a unit step function, fit) = 5(0. The response 
corresponding to an arbitrary motion fit) can then be written using Duhamel's 
superposition integral as 



Using the above equation, the response to a harmonic motion, e ' a> can 

determined. Since only the periodic response is of interest, the limit of the above 

integral as t is considered. Using a change of variable and extending the 
lower l imi t to - 00 , the following relation can be obtained 

m-iuFiMe** ( 12 ) 

where F^(w) is the Fourier transform of ( t ) given by 

-+00 

V 6 (co) = I F 5 it)e~ iwt dt (13) 


The indicial response method has the following drawback. The step change in 
the displacement results in infinite velocity at the time at which the step change 
occurs. The indicial response contains a large spike at this time. Therefore, very 
small time steps must be used to ensure that the results are accurate during this 
period of rapid transients and the accurate evaluation of the Fourier integral over 
this time interval is difficult. In addition, the treatment of this infinite velocity by 
finite differences leads to non-physical transients that can cause errors in the 
final results. To avoid these difficulties, researchers 23 have replaced the step 
function by a smooth version which does not result in infinite velocities and 
spikes. Polynomial functions of time have been used in place of the standard step 
function to obtain the necessary smoothness at the beginning and end of the step. 
This allows time steps of normal size to be used in the calculations. 

For an arbitrary motion fit) and the corresponding response Fit), the Fourier 
transform of Equation (11) gives 2 **: 

iuT b iu) = Fiu)lfiic) (14) 

where fioi) and Fiio) are the Fourier transforms of fit) and Fit), respectively. It 

should be noted from Equation (12) that icoF 5 iw) gives the response to a harmonic 
motion. Since only the response to a harmonic motion is to be obtained, any 

arbitrary motion and the corresponding response can be used to obtain ia)F b ia j) 
from Equation (14). 

Thus, the time variation of both the motion and the response are required to 
calculate the response to harmonic motions. To reduce the time required for the 
transients to decay, the smooth step is often replaced by a pulse 25 . In the pulse 
motion, the blade returns to its original position after the duration of the pulse. 
This is in contrast to the step motion in which the blade position is different before 
and after the step. The pulse motion thus allows the flowfield to return to its 
steady undisturbed state after the transients created by the pulse have decayed. 
The unsteady calculations therefore need to be carried out only long enough to 


13 


ensure that the solution has reached its final state (the same as the initial state) 
within some specified tolerance. 

Several smooth pulse shapes have been investigated by researchers. A 
comparison 3 of results from three different pulse shapes shows that although the 
shape and size of the pulse determines the transient response, the ratio of the 
Fourier transforms of the response and the motion remains unchanged. Thus, 
the shape of the pulse is not of particular importance although care must be taken 
to ensure that the transform of the pulse motion does not become zero in the 
frequency range of interest. In addition, the duration of the pulse must be selected 
according to the range of frequencies of interest. Thus, the harmonic time period 
of interest must be smaller than the duration of the pulse for the results to be 
meaningful. This places a lower limit on the values of frequency for which 
calculations can be made using a pulse of given duration. The upper limit on the 
frequency is determined by the size of the time step; the upper limit is generally 
not relevant because the time step is normally quite small for reasons of 
numerical stability and accuracy. 

In the present calculations, the pulse function is selected as: 
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t ttftrax 1 


fit) - 0 


for 0 £ t < tftiax 

for t £ tmax ( 15 ) 


where t max is the duration of the pulse. The above choice makes fit) and fit) 
vanish at t= 0 and t~ t max in addition, higher derivatives also go to zero at 
t— t max . This ensures that there is a smooth transition to and from the 
undisturbed blade position. 

The Pulse Response method is used in conjunction with the Influence 
Coefficient method as follows. One blade in the cascade is given a transient 

motion of the form hit) = h 0 fit) or ait) = a 0 fit). The calculations start with the 
steady solution and unsteady response to the pulse in either motion, plunging or 
pitching, is calculated until the transient flowfield reaches the steady flowfield 
within a specified tolerance. The motion as well as the responses on all the blades 
are recorded and Fourier transforms of these are calculated numerically for the 

frequency of interest. Using these transforms, the influence coefficients iQk,o) are 
calculated from Equations (14); it is to be noted that the harmonic response, 

iwF^iio) , obtained from Equation (14) for this case, is simply the influence 
coefficient, Qk,0 • Equation (10) is then used to calculate the frequency domain 
unsteady aerodynamic coefficients Hkh > lha » lah ^d la a) f° r the interblade 
phase angle of interest. In this way, the coefficients can be determined for various 
values of reduced frequency by recalculating the Fourier transforms for the 
frequency of interest using the same time histories. 
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RESULTS AND DISCUSSION 


The results are presented in two parts. The results in the initial part serve to 
validate the Influence Coefficient method and the Pulse Response method by 
comparison with the Harmonic Oscillation method. In the remaining part of this 
section, results of flutter calculations for two examples are presented. 

Calculation of Aerodynamic Coefficients 

Figure 2 shows the unsteady lift coefficient due to plunging motion in a 
cascade comprised of eight flat plates. The cascade is unstaggered and has a 
unity gap-to-chord ratio. The Mach number of the flow at the inlet is Moo-0.5 and 

the reduced frequency of oscillation is fc c =l . 0 . A 41 *21 grid is used in the 
calculations with 41 points in the streamwise direction and 21 points in the 
direction of the stagger line in each of the 8 interblade passages; 2 1 points are 
located on each (upper and lower) surface of each airfoil; all grid points are 
uniformly distributed. 150 time steps/cycle of oscillation are used which 
corresponds to a nondimensional time step of 0.084 and the calculations are 
continued for three cycles of blade oscillation to allow non-periodic starting 
transients to decay. The amplitude of oscillation used is h 0 lc — 0.002 which is 
found to be small enough to yield a response which is linearly dependent on the 
amplitude of motion. 

The results obtained from the IC method and the HO method are represented 
as open symbols; the results from classical linear theory 4 are shown as dashed 
lines and denoted as 'linear [4]'. A comparison of the various results reveals the 
following. The HO method and the IC method yield virtually identical results. 
This is simply a verification of the superposition principle and the associated 
calculations. It implicitly confirms that the unsteady problem is linear for the 
amplitude of motion used in the calculations. Similar comparisons made with 
other airfoils and other cascade geometries have demonstrated the same level of 
agreement seen here. The results of the present calculations also agree quite well 
with the results from classical linear theory. For flat plate airfoils and small 
amplitude oscillations, the present calculations are expected to yield the same 
results as obtained from classical linear theory. It is to be noted that classical 
linear theory uses a semi- analytical formulation and the solutions show singular 
behavior near the leading edge. This feature is not completely captured by the 
relatively coarse 41*21 grid used in the present calculations. This may account for 
some of the differences observed between the results. Figure 3 shows the moment 
coefficient for the plunging motion. Figures 4 and 5 show the lift and moment 
coefficients due to pitching about midchord; the amplitude of pitching oscillations 

used in the calculations is a o = 0.2° and 210 time steps/cycle of oscillation are used. 
In all the above results, the agreement between the HO method and the IC method 
is extremely good. The agreement between the present results and those from 
classical linear theory ranges from fair to good. 
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Figure 6 shows the time-variation of the displacement and the resulting lift on 
all the blades of a 8-blade cascade of flat plates for use in the PR method. The grid 
and the time step used are the same as those used in the previous calculations 
with the HO method and the IC method. The motion of blade 5 is Bhown as a 
function of time; the other blades in the cascade remain stationary. The 
rpgYimnm displacement during the pulse is hjc = 0.002 and the duration of the 

pulse is t max =20tt. The lift on all blades is shown separately. As can be seen, the 
calculations are continued until all responses have returned to the initial 
undisturbed values. For the present calculation, the initial and final response 
levels correspond to zero lift since the airfoils are flat plates. If the calculations 
are for a configuration which has a non -zero value of steady state lift, then 
similar variations will be obtained for the unsteady component of the lift 
(instantaneous lift minus steady lift). The time histories are Fourier transformed 
and combined according to Equation (14) to obtain the influence coefficients 
( QN-k,0 ) which are combined according to Equation (10) to get the harmonic 
coefficients. 

Figure 7 shows the lift due to plunging obtained by applying the PR method 
and the IC method to the time histories of Figure 6. The solid lines are results 
from the PR+IC method and the open symbols represent results from the HO 
method. Results from classical linear theory are also shown for comparison. The 
agreement between the PR+IC method and the HO method is very good. This 
substantiates the validity of the PR method used in combination with the IC 
method. It also implicitly confirms the linearity of the unsteady problem which 
permits the use of the Duhamel superposition integral and the influence 
coefficient approach. The agreement between the present calculations and the 
results from classical linear theory is also quite good except at values of reduced 
frequency near the point of acoustic resonance. The comments made earlier 
regarding the singularity near the leading edge and the limitations of the grid 
used in the present study can be applied to account for some of the observed 
discrepancy. Figures 8-10 show the variation of the remaining coefficients for the 
same cascade and flow condition. As before, the agreement between the PR+IC 
method and the HO method is very good and the agreement between the present 
results and those from classical linear theory varies from fair to good. 

All the computations described here were performed on a CRAY X-MP 
computer. The calculations performed using the HO method required about 160 

CPU seconds for the case of pitching motion with fe c =1.0 and <7=45°. 
Approximately the same amount of time was required for the corresponding 
calculations using the IC method with fc c =1.0 which gave results for 

o=0°, 45°, 90°,..., 315°. The calculation of the time histories in the PR method 
required about 230 CPU seconds; the additional time required for calculating 
Fourier transforms was quite small and was not recorded. 
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Calculation of Flutter Boundaries 

Flutter results are presented for two examples. The first example is a single- 
degree-of-freedom (pitching) system that has been previously considered in 

References 2 and 12. The cascade has nine blades, a stagger angle of 6=45°, a gap- 
to-chord ratio of g/c= 1.0 and the elastic axis is located at the leading edge 
(a/t=-1.0). The stability of this system can be inferred entirely from the imaginary 
part of the moment coefficient, Im{Z a }. Flutter occurs if lm{l a « } ^ 0 in which 
case the oscillating cascade extracts energy from the fluid stream. Figure 11 
shows the moment coefficient (real and imaginary parts) at the different values of 
interblade phase angle that may arise at flutter (Equation (5)). The results are for 
an inlet Mach number of Moo -0.5 and a reduced frequency of k c =0.222. 
Calculations have been performed on three different grids with 41*21, 81*41, and 
121*61 points in each interblade passage; the number of grid points on each airfoil 
surface (upper and lower) is 21, 41 and 61, respectively. The results from classical 
linear theory are included for comparison. 

It can be seen from Figure 11 that a uniform refinement of the grid spacing by 
a factor of 3 results only in small changes in the value of the moment coefficient. 
Some discrepancy is also observed between the present results and the results 
obtained from classical linear theory; this difference is not always reduced by grid 
refinement. This indicates that some other cause exists for the observed 
discrepancy. Efforts are currently underway to resolve this difficulty; no 
explanation is offered at present. It should be noted that observed difference 
cannot be attributed to the use of the PR method or the IC method since these have 
been validated by comparison with the HO method. 

For the calculations performed using the 41*21 grid, a reduced frequency of 

kc= 0.222 results in neutral stability (lm{/ aa } = 0) of the 0=320° aeroelastic mode 
with all other modes being stable (Im{/ a « } < 0); this result was obtained after 
calculating the coefficient la a ®t different values of reduced frequency and 
interblade phase angle. Thus, the reduced frequency at flutter is k c f= 0.222 and the 

phase angle at flutter is af- 320°. The corresponding results from classical linear 
theory are fc c /=0.254 and <Tf= 320°. Similar calculations have been performed for 
other values of the inlet Mach number between 0.2 and 0.8 using a 81*41 grid in 
each interblade passage. Figure 12 shows the variation of flutter reduced 
frequency with Mach number. It can be seen that the reduced frequency at flutter 
increases with the inlet Mach number with a levelling-off at the higher values of 
Mach number. In addition to the flat plate airfoil, results are also shown for a 
double-circular-arc airfoil with thickness-to-chord ratio of 0.05 (5% DCA airfoil). 
These calculations have also been performed with a 81*41 grid in each interblade 
passage. In the present example, the effect of thickness is to increase the flutter 
reduced frequency at all the values of inlet Mach number. The interblade phase 

angle at flutter was found to be <rf= 320° for all the results shown. 
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The second example for which flutter calculations have been performed 
consists of a 5-blade cascade in which the typical section has two degrees-of- 
freedom. This example has been previously considered in Reference 2; the 
geometric and structural parameters in this example are representative of the 

SR5 propfan 13 . The cascade stagger angle is ft=10.7° and the gap-to-chord ratio is 
g/c=l .85. The airfoil section is from a NACA 16 series with a thickness-to-chord 
ratio of 0.03 and a design lift coefficient ,of 0.3; this airfoil, which is taken from the 
three-quarter span location of the SR6 propfan, is hereafter referred to as the SR5 
airfoil. The structural model for each blade is a two-degrees-of-freedom typical 
section with elastic axis at the leading edge (a^—1.0). The mass ratio is ^"1 15, the 
radius of gyration is r a “1.076, the offset between elastic axis and center of mass is 
x a =0.964 and the ratio of uncoupled natural frequencies in bending and torsion is 
(Oh/Wa =0.567. 

The following procedure has been used in the calculations. For each value of 
interblade phase angle obtained from Equation (6), the reduced frequency is varied 

until one of the two eigenvalues of Equation (7) displays neutral stability (Jf=0 I) 

while the other eigenvalue displays stability (jti<0). The reduced frequency thus 
obtained (k c f) is shown in Figure 13 for a flat plate airfoil. Calculations have been 
done using three different grids with 41*38, 61*57, and 81*75 points in each of the 
five interblade passages; the number of grid points on each airfoil surface (upper 
and lower) is 21, 31 and 41, respectively. Results from classical linear theory are 
included for comparison. It can be seen that reducing the grid spacing by a factor 
of 2 does not result in much change in the reduced frequency at flutter. Figure 14 
shows the flutter reduced velocity for the same case. It can be seen that the lowest 

flutter reduced velocity occurs at <7=288° and this is identified as the critical phase 
angle. The corresponding flutter reduced velocity is V* f m 7.82 for the 41*38 grid, 
V*f “7.71 for the 61*57 grid and V*f * 7.63 for the 81*75 grid; the result from 
classical linear theory is V*f m 7.01. 

Similar calculations have been performed for values of inlet Mach number 
between 0.2 and 0.8; a 61*57 grid is used. Calculations are done for a flat plate 
airfoil and the SR6 airfoil. The results are presented in Figures 16 and 16 as 
variations of flutter reduced velocity (V* f) and flutter frequency ratio {uifliOa) with 
Mach number. Figure 16 shows that the flutter reduced velocity decreases with 
increasing Mach number for both the flat plate airfoil and the SR6 airfoil; the 
variation is seen to be almost linear. Figure 16 shows that the flutter frequency 
ratio also decreases with increasing Mach number. It is noted that there is only a 
marginal difference between the results for the flat plate and the SR5 airfoil. 
However, it may be recalled that the SR6 airfoil has a small camber angle, a small 
thickness-to-chord ratio of 0.03 and the cascade gap-to-chord ratio is large 
(g/c=l. 85). Therefore, it may be concluded that there is not much steady flow 
deflection for this case and consequently the flutter results do not show much 
difference between the flat plate airfoil and the SR5 airfoil. 
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CONCLUDING REMARKS 


Two methods of calculating (linear) frequency domain unsteady aerodynamic 
coefficients from a time-marching full-potential solver have been developed and 
verified. The first method, the Influence Coefficient method, allows coefficients 
for different interblade phase angles to be calculated simultaneously. The second 
method, the Pulse Response method, allows coefficients for several oscillation 
frequencies to be calculated from the same transient response. When these two 
methods are combined, the aerodynamic coefficients for several combinations of 
phase angles and frequencies can be calculated at approximately the 
computational cost required for the calculation of a single phase angle and single 
frequency using the Harmonic Oscillation method. These methods have been 
verified individually and in combination by comparison with the Harmonic 
Oscillation method. For small amplitudes of motion, the unsteady flow problem is 
linear and therefore these methods give accurate results, as expected. 

Flutter calculations have been performed for two examples. The first example 
has only pitching degree of freedom. Calculations have been done over a range of 
subsonic Mach numbers using a flat plate airfoil and a 5% thick double -circular- 
arc airfoil. For both airfoils, the flutter reduced frequency is seen to increase with 
Mach number and the effect of airfoil thickness is to increase the flutter reduced 
frequency. The results obtained from the present calculations with a flat plate 
airfoil show some differences when compared with the results from the classical 
linear theory. These differences are not eliminated by grid refinement indicating 
that some other cause contributes to this discrepancy. However, there is no 
indication that the problem lies in the use of either the Influence Coefficient 
method or Pulse Response method. 

The second example for which flutter calculations have been performed is a 
two-degree-of-freedom system with geometric and structural parameters 
representative of the SR5 propfan. Once again, calculations have been performed 
over a subsonic range of Mach numbers. The flutter reduced velocity and the 
flutter frequency ratio are seen to decrease continuously with increasing Mach 
number; the decrease in flutter reduced velocity is almost linear with Mach 
number. The difference between the results for the flat plate and the SR5 airfoil 
are negligible. It is inferred that the combination of large gap-to-chord ratio, 
small thickness-to-chord ratio and small camber angle results in very little mean 
flow deflection and consequently there is very little change in the unsteady 
aerodynamic behavior due to the addition of loading. 

The present approach allows a unified analysis capability in which both time 
domain and frequency domain flutter calculations can be performed using the 
same time-marching algorithm. This will allow a direct comparison of results 
from linear and nonlinear flutter analyses without concern for differences in CFD 
algorithms, grids and other purely numerical factors. The computational 
efficiency that is provided by these methods will allow flutter calculations to be 
done more routinely than was previously possible with the Harmonic Oscillation 
method. The methods implemented here on a full-potential solver can just as 
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easily be applied to other time-marching analyses, such as those based on the 
Euler equations. Although the present work has been restricted to subsonic Mach 
numbers, transonic flow calculations can also be performed using the same 

methods. 


20 



REFERENCES 


1 Mehmed, 0., Kaza, K. R. V., Lubomski, J. F., and Kielb, R. E., "Bending- 
Torsion Flutter of a Highly Swept Advanced Turboprop," NASA TM 82975, 1982. 

2 Bakhle, M. A., Reddy, T. S. R., $eith, T. G., Jr., "Time Domain Flutter 
Analysis of Cascades Using a Full-Potential Solver," AIAA Paper 90-0984, Apr. 
1990. To be published in AIAA Journal, July 1991. 

3 Williams, M. H. and Ku, C. C., "Three Dimensional Full Potential 
Aerodynamic Method for the Aeroelastic Modeling of Propfans," AIAA Paper 90- 
1120, Apr. 1990. 

4 Smith, S. N., "Discrete Frequency Sound Generation in Axial Flow 
Turbomachines," British Aeronautical Research Council, London, ARC R&M 
No. 3709, 1971. 

3 Lane, F, "Supersonic Flow Past an Oscillating Cascade with Supersonic 
Leading-Edge Locus," Journal of the Aeronautical Sciences, Vol. 24, Jan. 1957, 
pp. 65-66. 

6 Whitehead, D. S. and Grant, R. J., "Force and Moment Coefficients for High 
Deflection Cascades," Proc. 2nd Inti. Symp. on Aeroelasticity in Turbomachines, 
(ed. P. Suter), Juris-Verlag Zurich, 1981, pp. 85-127. 

7 Verdon, J. M. and Caspar, J. R., "Development of a Linear Aerodynamic 
Analysis for Finite -Deflection Subsonic Cascades," AIAA Journal, Vol. 20, No. 9, 
Sep. 1982, pp. 1259-1267. 

8 Hall, K. C. and Crawley, E. F., "Calculation of Unsteady Flows in 
Turbomachinery Using the Linearized Euler Equations," AIAA Journal, Vol. 27, 
No. 6, June 1989, pp. 777-787. 

9 Kao, Y. F., "A Two-Dimensional Unsteady Analysis for Transonic and 
Supersonic Cascade Flows," Ph.D. Thesis, School of Aeronautics and 
Astronautics, Purdue University, West Lafayette, Indiana, May 1989. 

10 Huff, D. L., "Numerical Analysis of Flow Through Oscillating Cascade 
Sections," AIAA Paper 89-0437, Jan. 1989. 

11 Reddy, T. S. R., Bakhle, M. A., and Huff, D. L., "Flutter Analysis of a 
Supersonic Cascade in Time Domain Using an Euler Solver," NASA TM in 
preparation, 1990. 

12 Bakhle, M. A., Keith, T. G. Jr., and Kaza, K. R. V., "Application of a Full- 
Potential Solver to Bending-Torsion Flutter in Cascades," AIAA Paper 89-1386, 
Apr. 1989. 


21 


13 Reddy, T. S. R., Srivastava, R., and Kaza, K. R. V., "The Effects of Rotational 
Flow, Viscosity, Thickness, and Shape on Transonic Flutter Dip Phenomena, 
AIAA Paper 88-2348, Apr. 1988. 

14 Bisplinghoff, R. L. and Ashley, H., "Principles of Aeroelasticity," John Wiley 
and Sons, Inc., New York, 1962. 

15 Shankar, V., Ide, H., Gorski, J., and Osher, S., "A Fast Time-Accurate 
Unsteady Full Potential Scheme," AIAA Paper 86-1512, Aug. 1985. 


16 Hedstrom, G. W., "Nonreflecting Boundary Conditions for Nonlinear 
Hyperbolic Systems," Journal of Computational Physics, Vol. 30, 1979, pp. 222-237. 

17 Lane, F., "System Mode Shapes in the Flutter of Compressor Blade Rows," 
Journal of the Aeronautical Sciences, Vol. 23, Jan. 1956, pp. 54-66. 


18 Bendiksen, O., and Friedmann, P., "Coupled Bending-Torsion Flutter in 
Cascades,” AIAA Journal, Vol. 18, No. 2, Feb. 1980, pp. 194-201. 


19 Nixon D Tzuoo, K. L., and Ayoub, A., "Rapid Computation of Unsteady 
Transonic Cascade Flows," AIAA Journal, Vol. 25, No. 5, May 1987, pp. 760-762. 


20 Buffum, D. H., and Fleeter, S., "Aerodynamics of a Linear Oscillating 
Cascade," NASA TM 103250, Aug. 1990. 

21 Ballhaus, W. F. and Goorjian, P. M., "Computation of Unsteady Transonic 
Flows by the Indicial Method," AIAA Journal, Vol. 16, No. 2, Feb. 1978, pp. 117- 
124. 

22 Kerlick, G. D. and Nixon, D., "A High-Frequency TransonicSmall 
Disturbance Code for Unsteady Flows in a Cascade," AIAA Paper 82-0955, June 
1982. 


23 Mvers M. R. and Ruo, S. Y., "Calculation of Unsteady Aerodynamic 
Coefficients Using Transonic Time Domain Methods," AIAA Paper 83-0885, Apr. 
1983. 


24 Davies D. E. and Salmond, D. J., "Indicial Approach to Harmonic 
Perturbations in Transonic Flow," AZAA Journal, Vol. 18, No. 8, Aug. 1980, pp. 
1012-1014. 


25 Seidel, D. A., Bennett, R. M., and Whitlow, W., Jr "An Exploratory Study of 
Finite-Difference Grids for Transonic Unsteady Aerodynamics, AIAA Paper 83- 
0503, Jan. 1983. 


22 



Figure la: Typical section blade model with two degrees-of-freedom. 



periodic 


Figure lb: Cascade geometry and schematic of H-type grid. 
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Figure 2: Variation of lift due to plunging with phase angle. 
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Figure 3: Variation of moment due to plunging with phase angle. 
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Figure 4: Variation of lift due to pitching with phase angle. 
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Figure 5: Variation of moment due to pitching with phase angle. 
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Figure 6: Time history of displacement and lift in Pulse Response method. 
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Figure 7: Variation of lift due to plunging with reduced frequency. 
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Figure 8: Variation of moment due to plunging with reduced frequency. 
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Figure 9: Variation of lift due to pitching with reduced frequency. 
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Figure 11: Variation of moment due to pitching with interblade phase angle, 
example 1. 
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Figure 12: Variation of flutter reduced frequency with Mach number, example 1. 
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Figure 13: Variation of flutter reduced frequency with interblade phase angle, 
example 2. 
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Figure 14: Variation of flutter reduced velocity with interblade phase angle, 
example 2. 
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Figure 15: Variation of flutter reduced velocity with Mach number, example 2. 
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Figure 16: Variation of flutter frequency ratio with Mach number, example 2. 
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